This section is intended to introduce various aspects of the art, which may be associated with exemplary embodiments of the present techniques. This discussion is believed to assist in providing a framework to facilitate a better understanding of particular aspects of the present techniques. Accordingly, it should be understood that this section should be read in this light, and not necessarily as admissions of prior art.
Many parameter estimation, inversion and imaging methods compute a forward simulation applying a sequence of steps extrapolating forward in time from an initial state. The associated inverse or imaging method applies an adjoint operator stepping backwards in time from a final time state. The method may then generate information about the system by cross correlating the forward simulation with the adjoint simulation at the same time steps. The information generated may then be used to improve the simulation parameters used to fit the available data, as well for other purposes. As an example, the information may be used to generate a gradient or a Hessian relating changes in an objective function to changes in the simulation parameters. In another example, the information may be used to create an image as done in Reverse Time Depth Migration (RTM).
These techniques require the forward-in-time simulation to be accessed in reverse time order for the cross-correlation computations. However, the size of the information needed to save the forward simulation for all time steps often exceeds the available memory storage.
For example, pre-stack Reverse Time depth Migration (RTM) is commonly used to image subsurface earth structures when doing hydrocarbon exploration with seismic data. During RTM, the receiver wavefield recorded in a seismic experiment is propagated backwards in time and cross-correlated with the forward-in-time simulation of the source wavefield. The means that the forward-in-time source simulation must be accessed in reversed time order. The straightforward approach would be to compute and save all time steps of the forward-propagated source wavefield at all subsurface locations. This can be impractical due to the huge data volumes involved. As an example, for a standard test project, the Society of Exploration Geophysicists (SEG) Advanced Modeling (SEAM) project, this could mean keeping 40,000 time steps of an 84-gigabyte subsurface volume for a total of about 3.36 petabytes of data. Not only is the amount of storage required problematic, but the access times required to move that amount of data may also be problematic.
To decrease the storage needed, a forward simulation wavefield state, termed herein a full state checkpoint, may be defined as including all the information needed to either perform a cross-correlation or to restart a forward simulation from a given time step. Full state checkpoints can then be saved at a smaller number of carefully chosen time steps so that the forward simulation can be restarted from the saved checkpoints and propagated to desired time steps. (The term checkpoint may also be used to refer to the aforementioned carefully chosen time steps, as contrasted with data saved at such time steps.) Thus, the forward simulation can be recomputed at time steps in reverse time order as needed starting from the saved checkpoints. The checkpoint memory is reused for new checkpoints whenever it is no longer needed. The trade-off is performing more computation to minimize the input/output and storage requirements. This technique may be useful whenever the data sizes involved are much larger than the memory storage available.
Griewank and Walther proposed an “optimal” way to choose which time steps should be checkpointed given that a specified amount of memory is available in which to save the checkpointed forward simulation wavefield states. See, e.g., Griewank, A., “Achieving logarithmic growth of temporal and spatial complexity in reverse automatic differentiation,” 1 Optimization Methods and Software 35-54 (1992); Griewank, A., Evaluating Derivatives: Principles and Techniques of Algorithmic Differentiation, Society for Industrial and Applied Mathematics (Philadelphia, Pa., 2000); Griewank, A. and A. Walther, “Algorithm 799: An implementation of checkpointing for the reverse or adjoint mode of computational differentiation,” 26 ACM Transactions on Mathematical Software 19-45 (2000).
In the SEAM example, it may be assumed that there are 40,000 time steps in the forward simulation. Using storage or memory to save 100-source wavefield state buffer checkpoints at any one time, the forward simulation can be accessed in reverse order using 114,747 forward simulation time steps, using the Griewank checkpointing scheme. Thus, a factor of roughly 400 in storage space is saved at an expense of doing 2.87 times as much computation for the forward simulation. For the RTM application, re-computing time steps for the forward simulation from checkpoints can often be faster than doing the input/output required to access the stored forward simulation time steps, even if sufficient disk space is available to store all of the forward simulation time steps in the case of slow disk and fast computation.
Symes applied Griewank optimal checkpointing as an effective implementation strategy for reverse time migration. Symes, W. W., 2007, Reverse time migration with optimal checkpointing, 72 Geophysics (No. 5), P.SM213-SM221. The checkpointing method may be particularly important for RTM that includes the physics of wavefield attenuation (e.g., imaging using both P-waves and S-waves), as the time reversal strategy of forward simulation using wavefield checkpointing is always stable. In comparison, an alternative RTM implementation strategy of saving boundary values and the final wavefield state and doing reverse time extrapolation of the source wavefield simulation can be unstable. Further, the instability in this technique makes it an inappropriate application strategy when attenuation is included as part of the physics in the forward simulation.
The checkpointing method for doing time reversal of a forward simulation has much more general application than just for RTM. Applications are quite general and associated with any time-stepping forward simulator. These include reservoir simulators, fluid flow, heat transfer, geologic basin modeling, as well as seismic full wavefield inversion (FWI). See, e.g., Tarantola, A., 1984, Inversion of seismic reflection data in the acoustic approximation: 49 Geophysics 1259-1266; Tarantola, A., Inverse Problem Theory: Method for Data Fitting and Model Parameter Estimation, Elsevier 125-133, 144-258 (1987); Plessix, R. E., “A review of the adjoint-state method for computing the gradient of a functional with geophysical applications,” 167 Geophysical Journal International 495-503 (2006). Krebs, J. R., J. E. Anderson, D. Hinkley, R. Neelamani, S. Lee, A. Baumstein, M. D. Lacasse, “Fast Full-Wavefield Seismic Inversion Using Encoded Sources,” 74 Geophysics P.WCC177-WCC188 (2009). The techniques may also apply to methods of automatic differentiation. See, e.g., Griewank, A., Juedes, D., and Srinivasan, “ADOL-C, a package for the automatic differentiation of algorithms written in C/C++,” Preprint MCS-180-1190, Mathematics and Computer Science Division, Argonne National Laboratory, Argonne, Ill. (1990); Griewank, A., Evaluating Derivatives: Principles and Techniques of Algorithmic Differentiation, Society for Industrial and Applied Mathematics (Philadelphia, Pa., 2000).